

*************************************************************************
*                                                                       *
*  Rotina para solucao de sistemas de equacoes lineares nao-simetricos  *
*                                                                       *
*  pelo metodo iterativo Generalized Minimum Residual - GMRES.          *
*                                                                       *
*  Esta implementacao se baseia no artigo de SHAKIB,F. , HUGES,J.R. ,   *
*                                                                       *
*  JOHAN,Z."A Multi-Element Group Preconditioned GMRES Algorithm for    *
*                                                                       *
*  nonsymetric Systems Arrising in Element Analysis" , Computer Methods *
*                                                                       *
*  in Applied Mechanics and Engineering, No 75 ,pags. 415-456, North    *
*                                                                       *
*  Holland, 1989.                                                       *
*                                                                       *
*************************************************************************



        SUBROUTINE GMRES(A,B,XIS,XISA,U,EBAR,EEBAR,HBAR,YIP,CE,
     &                 ES,ETOL,N,NB,K,LMAX,ICONV,ITER,IRP,IWR)


      IMPLICIT DOUBLE PRECISION (A-H,O-Z)



      DIMENSION A(N,N),B(N)
      DIMENSION XIS(N),XISA(N),U(N,NB+1),HBAR(NB+1,NB+1)
      DIMENSION YIP(NB+1),EBAR(NB+1),EEBAR(NB+1),CE(NB+1),ES(NB+1)
        
        dimension baux(1000)

C      COMMON /A/ D(2,2),XI(6,3),W(6,3),IDUP(500),INC(500,2),C(500),
C     *           S(500,3),ISYM(500),X(500),Y(500),IFIP(1000),    
C     *           A(1000,1000),P(1000),B(1000)C
C
C
C      COMMON /GMRS/ XIS(1000),U(1000,91),EBAR(91),HBAR(91,91),
C     *          YIP(91),CE(91),ES(91),           XISA(1000)


C...Descricao das Variaveis

C   N      => Numero de incognitas do sistema.
C   A(N,N) => Matriz de coeficientes do sistema.
C   B(N)   => Vetor independente do sistema.
C   ETOL   => Coeficiente arbitrado para estipular a tolerancia.
C   LMAX   => Numero maximo de iteracoes do GMRES.
C   EPSON  => Tolerancia.
C   XIS(N)   => Vetor que armazena a solucao aproximada.
C   U(N,K) => Matriz que armazena a base ortogonalizada para Z.
C   EBAR(K+1) => Vetor que vai armazenar o residuo na posicao k+1.
C   HBAR(K+1,K+1) => Matriz do sistema que resolvido vai gerar os coef.
C                    das bases ortogonalizadas.
C   YIP(K) => Vetor que acumula os coeficientes das bases ortogonais.




            r=0.
            do i=1,N
                    r = r + b(i)**2
            enddo
            write(*,*)'Antes solve gmres |b|= ',r


C...Inicializacao

       K1= K+1

C......Calcula a norma de B e zera X.
       XNORMX= 0.D0
       XNORMB= 0.D0
       DO 10 IAUX=1,N
          XNORMB= XNORMB + B(IAUX)*B(IAUX)
          XIS(IAUX)= 0.D0
          XISa(IAUX)= 0.D0
 10    CONTINUE
       XNORMB= DSQRT( XNORMB )

       EPSON= ETOL*XNORMB


c       WRITE (6,1000) EPSON
 1000  FORMAT(/,10X,'Norma do residuo admissivel = ',E12.6,/)

c       write(6,1010) XNORMB
 1010  format(/,5x,'Norma inicial do residuo = ',E12.6,/)


C...Ciclos do GMRES

      DO 20 L=1,LMAX
         ll=l

c         !  write(*,1030)
         IF(IRP.GE.1) THEN
           write(IWR,1030) 
 1030      format(/,10X,'Residuo nas iteracoes:'/
     *            /,17x,'Relativo',9x,'Absoluto',12x,
     *            '!!x!!',14x,'dr/dx',/)
         END IF

C...Calculo do primeiro vetor da base ortogonalizada de Z.

C......Zera U(N,K+1)
       DO 25 JAUX=1,K1
          DO 25 IAUX=1,N
             U(IAUX,JAUX)= 0.D0
 25     CONTINUE



C.......Calcula u1= b - A.X

        DO 30 JAUX=1,N
          DO 30 IAUX=1,N

              U(IAUX,1)= U(IAUX,1) - A(IAUX,JAUX)*XIS(JAUX)
 30     CONTINUE

        DO 35 IAUX=1,N

          U(IAUX,1)= U(IAUX,1) + B(IAUX)
 35     CONTINUE

C.......Calcula a norma de U e acumula em EBAR(1).

        EBAR(1)= 0.D0
        DO 50 IAUX=1,N
          EBAR(1)= EBAR(1) + U(IAUX,1)*U(IAUX,1)
 50     CONTINUE
        EBAR(1)= DSQRT( EBAR(1) )

        DO 60 IAUX=1,N
           U(IAUX,1)= U(IAUX,1) / EBAR(1)
 60     CONTINUE


C...Iteracao do GMRES.

        DO 70 I=1,K

           II=I
           I1=I+1

C.........Calcula U(i+1)= A.U(i)

          DO 80 JAUX=1,N
            DO 90 IAUX=1,N

              U(IAUX,I1)= U(IAUX,I1) + A(IAUX,JAUX)*U(JAUX,I)
 90         CONTINUE
 80       CONTINUE



C...Iteracao do algoritmo de ortogonalizacao de Gram-Schmidt modificado.

          DO 100 J=1,I

            HBAR(J,I)=0.D0
            DO 110 IAUX=1,N
               HBAR(J,I)= HBAR(J,I) + U(IAUX,I1)*U(IAUX,J)
 110        CONTINUE

            DO 120 IAUX=1,N
               U(IAUX,I1)= U(IAUX,I1) - HBAR(J,I)*U(IAUX,J)
 120        CONTINUE

 100      CONTINUE


C.........Calcula a norma de U(i+1).

          HBAR(I1,I)=0.D0
          DO 130 IAUX=1,N
            HBAR(I1,I)= HBAR(I1,I) + U(IAUX,I1)*U(IAUX,I1)
 130      CONTINUE

          HBAR(I1,I)= DSQRT( HBAR(I1,I) )

C.........Normaliza o vetor U(i+1).

          DO 140 IAUX=1,N
             U(IAUX,I1)= U(IAUX,I1)/HBAR(I1,I)
 140      CONTINUE


C.........Algoritmo Q-R

          IF(I.GT.1) THEN
            DO 150 J=1,I-1
              HBJI=   CE(J)*HBAR(J,I) + ES(J)*HBAR(J+1,I)
              HBJ1I= -ES(J)*HBAR(J,I) + CE(J)*HBAR(J+1,I)
              HBAR(J,I)= HBJI
              HBAR(J+1,I)= HBJ1I
 150        CONTINUE
          END IF


          R= HBAR(I,I)*HBAR(I,I) + HBAR(I1,I)*HBAR(I1,I)
          R= DSQRT(R)

          CE(I)= HBAR(I,I)/R
          ES(I)= HBAR(I1,I)/R

          HBAR(I,I)= R

          HBAR(I1,I)= 0.D0

          EBAR(I1)= -ES(I)*EBAR(I)
          EBAR(I)=   CE(I)*EBAR(I)

C........ Fim do algoritmo Q-R.

C...Teste de convergencia
          XNORMR= DABS(EBAR(I1))
          XNORXA= XNORMX
          inda=0

          if(i.eq.k) inda=1
          CALL CALCX(XIS,XISA,U,YIP,HBAR,EBAR,EEBAR,inda,I,N,NB,XNORMX)

C          CALL CALCX(inda,I,N,XNORMX)
           DRDX= (XNORMX-XNORXA)/(XNORMR-XNORRA)
          iter= (L-1)*K+I
          IF (ITER.NE.1) THEN
c            !  write(*,1060) iter,XNORMR/xnormb,XNORMR,XNORMX,DRDX
 1060       format(5x,i3,6X,e12.6,5X,e12.6,5X,E16.10,5X,E12.6)
            IF (IRP.GE.1) THEN
              write(IWR,1060) iter,XNORMR/xnormb,XNORMR,XNORMX,DRDX
            END IF
          ELSE
c            !  write(*,1061) iter,XNORMR/xnormb,XNORMR,XNORMX
 1061       format(5x,i3,6X,e12.6,5X,e12.6,5X,E16.10)
            IF (IRP.GE.1) THEN
              write(IWR,1061) iter,XNORMR/xnormb,XNORMR,XNORMX
            END IF
          END IF

          XNORRA= DABS(EBAR(I1))
          IF(DABS(EBAR(I1)).LE.EPSON) GOTO 1

 70     CONTINUE


C...Calculo dos Y por retro substituicao.

 1      I= II
c        CALL CALCX(I,N,XNORMX)
C...Teste de convergencia.

        IF(DABS(EBAR(I1)).LE.EPSON) then
          ICONV =1
          GOTO 2
        END IF

 20   CONTINUE

 2    L=LL

        if(iconv.eq.1)then
          iter= (L-1)*K+I

c          !  write(*,1070)I,L,iter
c          write(IWR,1070)I,L,iter
 1070     format(//,5x,'Convergencia atingida no passo:',i3,/,
     *              5X,'                      do ciclo:',i3,/,
     *              5X,'  ( ',I4,' iteracoes )',/)

        else
          !  write(*,1090)l
          write(IWR,1090)l
 1090     format(//,5x,'NAO foi atingida a convergencia em ',
     *              i3,' ciclos.')
        end if

C.......Calcula u1= b - A.X
        DO 195 IAUX=1,N
           U(IAUX,1)= 0.D0
 195    CONTINUE

        DO 200 JAUX=1,N
          DO 200 IAUX=1,N
              U(IAUX,1)= U(IAUX,1) - A(IAUX,JAUX)*XIS(JAUX)
 200     CONTINUE

        DO 210 IAUX=1,N
          U(IAUX,1)= U(IAUX,1) + B(IAUX)
 210     CONTINUE


C.......Calcula a norma de U e acumula em EBAR(1).

        RESID= 0.D0
        DO 220 IAUX=1,N
          RESID= RESID + U(IAUX,1)*U(IAUX,1)
 220    CONTINUE
        RESID= DSQRT( RESID )

        WRITE(IWR,1100) RESID
 1100   FORMAT(/,5X,'Norma do residuo calculado por  RESID = b - A.x',
     *        //,5x,'    Norma: ',e12.6,/)


C...ARMAZENA OS RESULTADOS NO VETOR B
        DO 190 IAUX=1,N
          B(IAUX)= XIS(IAUX)
 190    CONTINUE

        do i=1,N
            baux(i) = 0.
        enddo

        do i=1,N
            do j=1,N
                baux(i) = baux(i) + b(j)*A(i,j)
            enddo
        enddo


            r=0.
            do i=1,N
                    r = r + baux(i)**2
            enddo
            write(*,*)'depois solve gmres |A.x|= ',r


            r=0.
            do i=1,N
                    r = r + b(i)**2
            enddo
            write(*,*)'depois solve gmres |b|= ',r

            write(*,*)'iteracoes: ',iter,'equacoes: ',N

!         do i=1,N
!             write(*,*)baux(i)
!         enddo
        write(*,*)

      RETURN
      END
